#include <StVenKir.hpp>
#include <iostream>
#include <life/lifefilters/geotool.hpp>

namespace Life {

template<int Dim, int Order, template<uint16_type,uint16_type,uint16_type> class Entity>
typename StVenantKirchhoff<Dim,Order,Entity>::mesh_ptrtype
StVenantKirchhoff<Dim,Order,Entity>::createMesh( double meshSize )
{

    mesh_ptrtype mesh( new mesh_type );

    //Geometry
    GeoTool::Node x1(0,0);
    GeoTool::Node x2(4,0.1);
    GeoTool::Rectangle R( meshSize,"Structure",x1,x2);
    R.setMarker(_type="line",_name="Libre",_marker3=true);
    R.setMarker(_type="line",_name="Force",_marker1=true);
    R.setMarker(_type="line",_name="Fixe",_marker2=true,_marker4=true);
    R.setMarker(_type="surface",_name="OmegaStructure",_marker1=true);

this->application()->changeRepository( boost::format( "%1%/%2%/P%3%/h_%4%/" )
                                             % this->application()->about().appName()
                                             % convex_type::name()
                                             % Order
                                             % this->application()->vm()["hsize"].template as<double>()
                                             );
#if 0
    this->application()->changeRepository( boost::format( "%1%/%2%/P%3%/h_%4%/" )
                                           % this->application()->about().appName()
                                           % convex_type::name()
                                           % Order
                                           % this->application()->vm()["hsize"].template as<double>()
                                           );
#endif
    Gmsh gmsh;
    std::string fname = gmsh.generate( "domain", R.geoStr() );

    ImporterGmsh<mesh_type> import( fname );
    import.setVersion( "2.1" );
    mesh->accept( import );

    return mesh;

} // StVenantKirchhoff::createMesh

template<int Dim, int Order, template<uint16_type,uint16_type,uint16_type> class Entity>
void
StVenantKirchhoff<Dim, Order, Entity>::solve( sparse_matrix_ptrtype const& D,
                                   element_struct_type& u,
                                   vector_ptrtype const& F,
                                   bool is_sym )
{
    vector_ptrtype U( M_backend->newVector( u.functionSpace() ) );
    M_backend->solve( D, D, U, F, false );
    u = *U;
} // StVenantKirchhoff::solve

template<int Dim, int Order, template<uint16_type,uint16_type,uint16_type> class Entity>
void
StVenantKirchhoff<Dim, Order, Entity>::exportResults( element_struct_type& U, element_struct_type& V )
{
    M_exporter->step( 0 )->setMesh( U.functionSpace()->mesh() );
    M_exporter->step( 0 )->add( "u", U );
    M_exporter->save();
} // StVenantKirchhoff::export


template<int Dim, int Order, template<uint16_type,uint16_type,uint16_type> class Entity>
void
StVenantKirchhoff<Dim, Order, Entity>::exportResultsInstat( element_struct_type& U, element_struct_type& V,double time )
{
    M_exporter->step( M_cpt_export )->setMesh( U.functionSpace()->mesh() );
    M_exporter->step( M_cpt_export )->add( "u", U );
    M_exporter->save();

    ++M_cpt_export;

} // StVenantKirchhoff::export



} // Life

// instantiation
template class Life::StVenantKirchhoff<2,2,Life::Simplex>;
//template class Life::StVenantKirchhoff<3,2,Life::Simplex>;

